
############################### GSEA 富集分析 ##################################
library(msigdbr)
library(Seurat)
library(fgsea)
library(clusterProfiler) # GSEA富集/基因集读取
# BiocManager::install("clusterProfiler")
setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")
rm(list = ls())

#CD4_Treg####
# 基因集
geneSet_onco <- read.gmt("data/input/h.all.v2023.2.Hs.symbols.gmt")
# geneSet_onco <- read.gmt("data/input/c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt")#
# geneSet_onco <- read.gmt("data/input/c2.all.v2023.2.Hs.symbols.gmt")#出来的很杂乱
head(geneSet_onco)

# 导入差异分析后的数据，以便后续使用logFC进行基因排序
mydata <- read.csv("data/temp/CD4_Treg_DEG_results_to_GSEA.csv",row.names=1)
head(mydata)
mydata <- mydata[order(mydata$log2FoldChange,decreasing=TRUE),]#按log从大到小排序
mydata <- mydata[mydata$log2FoldChange!=Inf & mydata$log2FoldChange!=-Inf ,]#去除无穷大的，不然后面报错，Not all stats values are finite numbers
FCgenelist <- mydata$log2FoldChange
names(FCgenelist) <- rownames(mydata)
head(FCgenelist)

# 开始GSEA富集分析
GSEA_enrichment <- GSEA(FCgenelist,                 # 排序后的gene
                        TERM2GENE = geneSet_onco, # 基因集
                        pvalueCutoff = 0.05,      # P值阈值
                        minGSSize = 20,           # 最小基因数量
                        maxGSSize = 1000,         # 最大基因数量
                        eps = 0,                  # P值边界
                        pAdjustMethod = "BH")     # 校正P值的计算方法

result <- data.frame(GSEA_enrichment)
dim(GSEA_enrichment@result)
write.csv(result, "data/temp/CD4_Treg_GSEA_result_H.csv", row.names = FALSE)
result$ID

# 可视化
library(enrichplot) # 富集结果可视化

# 特定通路作图——单个通路
p1 <- gseaplot2(GSEA_enrichment, "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", color = "red3", pvalue_table = T) 
p1

# 特定通路绘图——多个通路
p2 <- gseaplot2(GSEA_enrichment, 
                result$ID, 
                #c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", "HALLMARK_TNFA_SIGNALING_VIA_NFKB"), 
          #color = c("red3", "blue4","darkgreen"), 
          pvalue_table = T) 
p2

# 展示富集到的通路，我们这里选择展示前15个
p3 <- dotplot(GSEA_enrichment, 
              showCategory = 15, 
              color = "p.adjust",
              title="CD4_treg") 
p3

# 将通路分为激活和抑制两个部分进行展示
library(ggplot2)     # 画图图
p4 <- dotplot(GSEA_enrichment, showCategory = 15, split = ".sign") + 
  facet_grid(~.sign) +
  theme(plot.title = element_text(size = 10, color = "black", hjust = 0.5),
        axis.title = element_text(size = 10,color = "black"), 
        axis.text = element_text(size = 10,color = "black"),
        axis.text.x = element_text(angle = 0, hjust = 1 ),
        legend.position = "right",
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))+
  ggtitle("CD4_Treg")
p4
pdf("./data/output/CD4_Treg_GSEA_H_气泡图_富集GSEA图.pdf",width = 12,height = 6)
# p1
p2
p3
p4
dev.off()


#没看懂代码
library(tidyverse)
rank_Signal2Noise <- read_tsv("./gsea_data/ranked_Signal2Noise_gene_list_lapa_versus_ctrl_1698913651364.tsv")
head(rank_Signal2Noise)
rank_Signal2Noise <- rank_Signal2Noise[,-4]
saveRDS(rank_Signal2Noise, "./gsea_data/rank_Signal2Noise.rds")


rank_Signal2Noise <- readRDS("./gsea_data/rank_Signal2Noise.rds")
head(rank_Signal2Noise)
# # A tibble: 6 × 3
#   NAME    TITLE                                                                            SCORE
#   <chr>   <chr>                                                                            <dbl>
# 1 SLC15A2 solute carrier family 15 member 2 [Source:HGNC Symbol;Acc:HGNC:10921]             4.46
# 2 MUC19   mucin 19, oligomeric [Source:HGNC Symbol;Acc:HGNC:14362]                          4.32
# 3 ELSPBP1 epididymal sperm binding protein 1 [Source:HGNC Symbol;Acc:HGNC:14417]            4.29
# 4 AGXT    alanine--glyoxylate aminotransferase [Source:HGNC Symbol;Acc:HGNC:341]            4.28
# 5 CYP4F8  cytochrome P450 family 4 subfamily F member 8 [Source:HGNC Symbol;Acc:HGNC:2648]  3.91
# 6 MMP13   matrix metallopeptidase 13 [Source:HGNC Symbol;Acc:H

#_NKT####
# 基因集
geneSet_onco <- read.gmt("data/input/h.all.v2023.2.Hs.symbols.gmt")
# geneSet_onco <- read.gmt("data/input/c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt")#
# geneSet_onco <- read.gmt("data/input/c2.all.v2023.2.Hs.symbols.gmt")#出来的很杂乱
head(geneSet_onco)

# 导入差异分析后的数据，以便后续使用logFC进行基因排序
mydata <- read.csv("data/output/NKT_DEG_results_to_KMplot.csv",row.names=1)
head(mydata)
mydata <- mydata[order(mydata$log2FoldChange,decreasing=TRUE),]#按log从大到小排序
mydata <- mydata[mydata$log2FoldChange!=Inf,]#去除无穷大的，不然后面报错，Not all stats values are finite numbers
FCgenelist <- mydata$log2FoldChange
names(FCgenelist) <- rownames(mydata)
head(FCgenelist)

# 开始GSEA富集分析
GSEA_enrichment <- GSEA(FCgenelist,                 # 排序后的gene
                        TERM2GENE = geneSet_onco, # 基因集
                        pvalueCutoff = 0.05,      # P值阈值
                        minGSSize = 20,           # 最小基因数量
                        maxGSSize = 1000,         # 最大基因数量
                        eps = 0,                  # P值边界
                        pAdjustMethod = "BH")     # 校正P值的计算方法

result <- data.frame(GSEA_enrichment)
dim(GSEA_enrichment@result)
write.csv(result, "data/output/NKT_GSEA_result_H.csv", row.names = FALSE)
result$ID

# 可视化
library(enrichplot) # 富集结果可视化

# 特定通路作图——单个通路
p1 <- gseaplot2(GSEA_enrichment, "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", color = "red3", pvalue_table = T) 
p1

# 特定通路绘图——多个通路
p2 <- gseaplot2(GSEA_enrichment, 
                result$ID, #c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", "HALLMARK_TNFA_SIGNALING_VIA_NFKB"), 
                #color = c("red3", "blue4","darkgreen"),
                pvalue_table = T) 
p2

# 展示富集到的通路，我们这里选择展示前15个
p3 <- dotplot(GSEA_enrichment, showCategory = 15, color = "p.adjust") 
p3

# 将通路分为激活和抑制两个部分进行展示
library(ggplot2)     # 画图图
p4 <- dotplot(GSEA_enrichment, showCategory = 10, split = ".sign") + facet_grid(~.sign) +
  theme(plot.title = element_text(size = 10, color = "black", hjust = 0.5),
        axis.title = element_text(size = 10,color = "black"), 
        axis.text = element_text(size = 10,color = "black"),
        axis.text.x = element_text(angle = 0, hjust = 1 ),
        legend.position = "right",
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))
p4
pdf("./data/output/NKT_GSEA_H_气泡图_富集GSEA图.pdf",width = 12,height = 6)
p2
p3
p4
dev.off()


#没看懂代码
library(tidyverse)
rank_Signal2Noise <- read_tsv("./gsea_data/ranked_Signal2Noise_gene_list_lapa_versus_ctrl_1698913651364.tsv")
head(rank_Signal2Noise)
rank_Signal2Noise <- rank_Signal2Noise[,-4]
saveRDS(rank_Signal2Noise, "./gsea_data/rank_Signal2Noise.rds")


rank_Signal2Noise <- readRDS("./gsea_data/rank_Signal2Noise.rds")
head(rank_Signal2Noise)
# # A tibble: 6 × 3
#   NAME    TITLE                                                                            SCORE
#   <chr>   <chr>                                                                            <dbl>
# 1 SLC15A2 solute carrier family 15 member 2 [Source:HGNC Symbol;Acc:HGNC:10921]             4.46
# 2 MUC19   mucin 19, oligomeric [Source:HGNC Symbol;Acc:HGNC:14362]                          4.32
# 3 ELSPBP1 epididymal sperm binding protein 1 [Source:HGNC Symbol;Acc:HGNC:14417]            4.29
# 4 AGXT    alanine--glyoxylate aminotransferase [Source:HGNC Symbol;Acc:HGNC:341]            4.28
# 5 CYP4F8  cytochrome P450 family 4 subfamily F member 8 [Source:HGNC Symbol;Acc:HGNC:2648]  3.91
# 6 MMP13   matrix metallopeptidase 13 [Source:HGNC Symbol;Acc:H


#_CD8_Tex####
# 基因集
geneSet_onco <- read.gmt("data/input/h.all.v2023.2.Hs.symbols.gmt")
# geneSet_onco <- read.gmt("data/input/c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt")#
# geneSet_onco <- read.gmt("data/input/c2.all.v2023.2.Hs.symbols.gmt")#出来的很杂乱
head(geneSet_onco)

# 导入差异分析后的数据，以便后续使用logFC进行基因排序
mydata <- read.csv("data/output/CD8_Tex_DEG_results_to_KMplot.csv",row.names=1)
head(mydata)
mydata <- mydata[order(mydata$log2FoldChange,decreasing=TRUE),]#按log从大到小排序
mydata <- mydata[mydata$log2FoldChange!=Inf,]#去除无穷大的，不然后面报错，Not all stats values are finite numbers
FCgenelist <- mydata$log2FoldChange
names(FCgenelist) <- rownames(mydata)
head(FCgenelist)

# 开始GSEA富集分析
GSEA_enrichment <- GSEA(FCgenelist,                 # 排序后的gene
                        TERM2GENE = geneSet_onco, # 基因集
                        pvalueCutoff = 0.05,      # P值阈值
                        minGSSize = 20,           # 最小基因数量
                        maxGSSize = 1000,         # 最大基因数量
                        eps = 0,                  # P值边界
                        pAdjustMethod = "BH")     # 校正P值的计算方法

result <- data.frame(GSEA_enrichment)
dim(GSEA_enrichment@result)
write.csv(result, "data/output/CD8_Tex_GSEA_result_H.csv", row.names = FALSE)
result$ID

# 可视化
library(enrichplot) # 富集结果可视化

# 特定通路作图——单个通路
p1 <- gseaplot2(GSEA_enrichment, "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", color = "red3", pvalue_table = T) 
p1

# 特定通路绘图——多个通路
p2 <- gseaplot2(GSEA_enrichment, 
                result$ID, #c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", "HALLMARK_TNFA_SIGNALING_VIA_NFKB"), 
                #color = c("red3", "blue4","darkgreen"),
                pvalue_table = T) 
p2

# 展示富集到的通路，我们这里选择展示前15个
p3 <- dotplot(GSEA_enrichment, showCategory = 15, color = "p.adjust") 
p3

# 将通路分为激活和抑制两个部分进行展示
library(ggplot2)     # 画图图
p4 <- dotplot(GSEA_enrichment, showCategory = 10, split = ".sign") + facet_grid(~.sign) +
  theme(plot.title = element_text(size = 10, color = "black", hjust = 0.5),
        axis.title = element_text(size = 10,color = "black"), 
        axis.text = element_text(size = 10,color = "black"),
        axis.text.x = element_text(angle = 0, hjust = 1 ),
        legend.position = "right",
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))
p4
pdf("./data/output/CD8_Tex_GSEA_H_气泡图_富集GSEA图.pdf",width = 12,height = 6)
p2
p3
p4
dev.off()


#没看懂代码
library(tidyverse)
rank_Signal2Noise <- read_tsv("./gsea_data/ranked_Signal2Noise_gene_list_lapa_versus_ctrl_1698913651364.tsv")
head(rank_Signal2Noise)
rank_Signal2Noise <- rank_Signal2Noise[,-4]
saveRDS(rank_Signal2Noise, "./gsea_data/rank_Signal2Noise.rds")


rank_Signal2Noise <- readRDS("./gsea_data/rank_Signal2Noise.rds")
head(rank_Signal2Noise)
# # A tibble: 6 × 3
#   NAME    TITLE                                                                            SCORE
#   <chr>   <chr>                                                                            <dbl>
# 1 SLC15A2 solute carrier family 15 member 2 [Source:HGNC Symbol;Acc:HGNC:10921]             4.46
# 2 MUC19   mucin 19, oligomeric [Source:HGNC Symbol;Acc:HGNC:14362]                          4.32
# 3 ELSPBP1 epididymal sperm binding protein 1 [Source:HGNC Symbol;Acc:HGNC:14417]            4.29
# 4 AGXT    alanine--glyoxylate aminotransferase [Source:HGNC Symbol;Acc:HGNC:341]            4.28
# 5 CYP4F8  cytochrome P450 family 4 subfamily F member 8 [Source:HGNC Symbol;Acc:HGNC:2648]  3.91
# 6 MMP13   matrix metallopeptidase 13 [Source:HGNC Symbol;Acc:H


#NKT####
# 基因集
geneSet_onco <- read.gmt("data/input/h.all.v2023.2.Hs.symbols.gmt")
# geneSet_onco <- read.gmt("data/input/c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt")#
# geneSet_onco <- read.gmt("data/input/c2.all.v2023.2.Hs.symbols.gmt")#出来的很杂乱
head(geneSet_onco)

# 导入差异分析后的数据，以便后续使用logFC进行基因排序
mydata <- read.csv("data/temp/NKT_DEG_results_to_GSEA.csv",row.names=1)
head(mydata)
mydata <- mydata[order(mydata$log2FoldChange,decreasing=TRUE),]#按log从大到小排序
mydata <- mydata[mydata$log2FoldChange!=Inf & mydata$log2FoldChange!=-Inf ,]#去除无穷大的，不然后面报错，Not all stats values are finite numbers
FCgenelist <- mydata$log2FoldChange
names(FCgenelist) <- rownames(mydata)
head(FCgenelist)

# 开始GSEA富集分析
GSEA_enrichment <- GSEA(FCgenelist,                 # 排序后的gene
                        TERM2GENE = geneSet_onco, # 基因集
                        pvalueCutoff = 0.05,      # P值阈值
                        minGSSize = 20,           # 最小基因数量
                        maxGSSize = 1000,         # 最大基因数量
                        eps = 0,                  # P值边界
                        pAdjustMethod = "BH")     # 校正P值的计算方法

result <- data.frame(GSEA_enrichment)
dim(GSEA_enrichment@result)
write.csv(result, "data/temp/NKT_GSEA_result_H.csv", row.names = FALSE)
result$ID

# 可视化
library(enrichplot) # 富集结果可视化

# 特定通路作图——单个通路
p1 <- gseaplot2(GSEA_enrichment, "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", color = "red3", pvalue_table = T) 
p1

# 特定通路绘图——多个通路
p2 <- gseaplot2(GSEA_enrichment, 
                result$ID, 
                #c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", "HALLMARK_TNFA_SIGNALING_VIA_NFKB"), 
                #color = c("red3", "blue4","darkgreen"), 
                pvalue_table = T) 
p2

# 展示富集到的通路，我们这里选择展示前15个
p3 <- dotplot(GSEA_enrichment, 
              showCategory = 30, 
              # color = "p.adjust",
              color = "GeneRatio",
              # orderBy = 'NES',
              # x = 'enrichmentScore',
              title="NKT") 
p3

# 将通路分为激活和抑制两个部分进行展示
library(ggplot2)     # 画图图
p4 <- dotplot(GSEA_enrichment, showCategory = 30, split = ".sign") + 
  facet_grid(~.sign) +
  theme(plot.title = element_text(size = 10, color = "black", hjust = 0.5),
        axis.title = element_text(size = 10,color = "black"), 
        axis.text = element_text(size = 10,color = "black"),
        axis.text.x = element_text(angle = 0, hjust = 1 ),
        legend.position = "right",
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 10))+
  ggtitle("NKT")
p4
pdf("./data/output/NKT_GSEA_H_气泡图_富集GSEA图.pdf",width = 12,height = 6)
# p1
p2
p3
p4
dev.off()


#没看懂代码
library(tidyverse)
rank_Signal2Noise <- read_tsv("./gsea_data/ranked_Signal2Noise_gene_list_lapa_versus_ctrl_1698913651364.tsv")
head(rank_Signal2Noise)
rank_Signal2Noise <- rank_Signal2Noise[,-4]
saveRDS(rank_Signal2Noise, "./gsea_data/rank_Signal2Noise.rds")


rank_Signal2Noise <- readRDS("./gsea_data/rank_Signal2Noise.rds")
head(rank_Signal2Noise)
# # A tibble: 6 × 3
#   NAME    TITLE                                                                            SCORE
#   <chr>   <chr>                                                                            <dbl>
# 1 SLC15A2 solute carrier family 15 member 2 [Source:HGNC Symbol;Acc:HGNC:10921]             4.46
# 2 MUC19   mucin 19, oligomeric [Source:HGNC Symbol;Acc:HGNC:14362]                          4.32
# 3 ELSPBP1 epididymal sperm binding protein 1 [Source:HGNC Symbol;Acc:HGNC:14417]            4.29
# 4 AGXT    alanine--glyoxylate aminotransferase [Source:HGNC Symbol;Acc:HGNC:341]            4.28
# 5 CYP4F8  cytochrome P450 family 4 subfamily F member 8 [Source:HGNC Symbol;Acc:HGNC:2648]  3.91
# 6 MMP13   matrix metallopeptidase 13 [Source:HGNC Symbol;Acc:H


# GSEA####
library(msigdbr)
library(Seurat)
library(fgsea)
library(clusterProfiler) # GSEA富集/基因集读取
# BiocManager::install("clusterProfiler")
setwd("C:/Users/ZFB/Desktop/单细胞生信/GSE212966")
rm(list = ls())

# 基因集
library(ggplot2)  
library(enrichplot) # 富集结果可视化
# devtools::install_github("junjunlab/GseaVis")
library(GseaVis)
geneSet_onco <- read.gmt("data/input/h.all.v2023.2.Hs.symbols.gmt")
# geneSet_onco <- read.gmt("data/input/c2.cp.kegg_legacy.v2023.2.Hs.symbols.gmt")#
# geneSet_onco <- read.gmt("data/input/c2.all.v2023.2.Hs.symbols.gmt")#出来的很杂乱
head(geneSet_onco)


##1. T sub_循环_GSEA####
sub=c("CD4_Tn", "CD4_Th", "CD4_Treg", "CD4_Tem", 
      "NKT", "CD8_Tc", "CD8_Te", "CD8_Tem", "CD8_Tex", "CD8_Trm")

set.seed(123)
for(sub_ in sub){
# sub_=sub[5]
  read_path <- paste0("./data/temp/",sub_,"_DEG_results_to_GSEA.csv")
  save_path_csv <- paste0("./data/temp/",sub_,"_GSEA_result_H.csv")
  save_path_polt <- paste0("./data/output/",sub_,"_GSEA_result_H_GSEA通路富集.pdf")
  read_path
  save_path_csv
  save_path_polt
  # 导入差异分析后的数据，以便后续使用logFC进行基因排序
  mydata <- read.csv(read_path ,row.names=1)
  head(mydata)
  mydata <- mydata[order(mydata$log2FoldChange,decreasing=TRUE),]#按log从大到小排序
  mydata <- mydata[mydata$log2FoldChange!=Inf & mydata$log2FoldChange!=-Inf ,]#去除无穷大的，不然后面报错，Not all stats values are finite numbers
  FCgenelist <- mydata$log2FoldChange
  names(FCgenelist) <- rownames(mydata)
  head(FCgenelist)
  # 开始GSEA富集分析
  GSEA_enrichment <- GSEA(FCgenelist,               # 排序后的gene
                          TERM2GENE = geneSet_onco, # 基因集
                          pvalueCutoff = 0.05,      # P值阈值
                          minGSSize = 20,           # 最小基因数量
                          maxGSSize = 1000,         # 最大基因数量
                          eps = 0,                  # P值边界
                          pAdjustMethod = "BH")     # 校正P值的计算方法,可以做P矫正BH法，这样去掉假阳性
  #遇到报错重启RStudio解决：GSEA analysis... Error in serialize(data, node$con) : error writing to connection
  
  result <- data.frame(GSEA_enrichment)
  dim(GSEA_enrichment@result)
  write.csv(result, save_path_csv, row.names = FALSE)
  
  result$ID
  GSEA_enrichment@result$Description<- sub("HALLMARK_", "", GSEA_enrichment@result$ID)
  
  # 可视化
  # 特定通路作图——单个通路
  # p1 <- gseaplot2(GSEA_enrichment, "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", color = "red3", pvalue_table = T) 
  # p1
  
  # 特定通路绘图——多个通路
  # p2 <- gseaplot2(GSEA_enrichment,
  #                 # result$ID,
  #                 c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", 
  #                   "HALLMARK_TNFA_SIGNALING_VIA_NFKB",
  #                   "HALLMARK_HYPOXIA"),
  #                 #color = c("red3", "blue4","darkgreen"),
  #                 pvalue_table = T)
  # p2
  #单独通路的GSEA图
  # mygene<- GSEA_enrichment@geneSets[["HALLMARK_HYPOXIA"]]
  #   # c("LOX","ANGPTL4","SDC4","CCN1","COL5A1","GPC1","SDC2","BGN",
  #   #                  "GPC1",
  #   #                  "GPC3",
  #   #                  "GPC4")
  # p21 <- gseaNb(object = GSEA_enrichment,
  #        geneSetID = 'HALLMARK_HYPOXIA',
  #        newGsea = T,
  #        # addPoint = F,
  #        # newHtCol = c("blue","white", "red")
  #        addGene = mygene
  #        )
  # P21
  
  # mygene<- GSEA_enrichment@geneSets[["HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"]]
  # p22 <- gseaNb(object = GSEA_enrichment,
  #               geneSetID = 'HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION',
  #               newGsea = T,
  #               # addPoint = F,
  #               # newHtCol = c("blue","white", "red")
  #               addGene = mygene
  # )
  # p22
  # 展示富集到的通路，我们这里选择展示前15个
  p3 <- dotplot(GSEA_enrichment, 
                showCategory = 50,
                # color = "p.adjust",
                color = "p.adjust",
                # orderBy = 'NES',
                x = "NES",
                # x = 'enrichmentScore',
                title=sub_) 
  # p3
  
  # 将通路分为激活和抑制两个部分进行展示
  # 画图
  p4 <- dotplot(GSEA_enrichment, showCategory = 30, split = ".sign") + 
    facet_grid(~.sign) +
    theme(plot.title = element_text(size = 10, color = "black", hjust = 0.5),
          axis.title = element_text(size = 10,color = "black"), 
          axis.text = element_text(size = 10,color = "black"),
          axis.text.x = element_text(angle = 0, hjust = 1 ),
          legend.position = "right",
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 10))+
    ggtitle(sub_)
  # p4
  
  #山脊图
  p5 <- ridgeplot(
    GSEA_enrichment,
    showCategory = 30,
    fill = "p.adjust",
    core_enrichment = TRUE,
    label_format = 30,
    orderBy = "NES",
    decreasing = FALSE)+
    ggtitle(sub_)
  # p5
  #气泡图-棒棒糖
  p6 <- dotplotGsea(data = GSEA_enrichment,
                    #topn = 30,
                    order.by = 'NES',
                    # line.type = 'dashed',
                    # pval=0.05, #defalut=0.05
                    # line.col = 'orange',
                    # add.seg = T
                    )
  # p6
  p7 <- dotplotGsea(data = GSEA_enrichment,
                    #topn = 30,
                    order.by = 'NES',
                    # line.type = 'dashed',
                    # pval=0.05, #defalut=0.05
                    # line.col = 'orange',
                    add.seg = T)
  # p7
  
  save_path_polt <- paste0("./data/output/",sub_,"_GSEA_result_H_GSEA通路富集8x8.pdf")
  pdf(save_path_polt,width = 8,height = 8)
  # p1
  # p2
  print(p3)
  print(p4)
  print(p5)
  print(p6)
  print(p7)
  #不用printpdf打开显示文件损坏
  dev.off()
  save_path_polt <- paste0("./data/output/",sub_,"_GSEA_result_H_GSEA通路富集12x8.pdf")
  pdf(save_path_polt,width = 12,height = 8)
  # p1
  # p2
  print(p3)
  print(p4)
  print(p5)
  print(p6)
  print(p7)
  dev.off()
}
#barplot和dotplot都只显示最显著的富集项，而用户可能想知道哪些基因与这些显著项有关。为了考虑基因可能属于多个注释类别的潜在生物学复杂性，并提供可用的数值变化信息，我们开发了cnetplot函数来提取复杂的关联。cnetplot将基因和生物学概念(例如GO terms或KEGG pathways)之间的联系描述为一个网络。GSEA结果也只支持核心富集基因的表达。
# mygene<- GSEA_enrichment@geneSets[["HALLMARK_HYPOXIA"]]
# cnetplot(GSEA_enrichment, foldChange=FCgenelist, circular = TRUE, colorEdge = TRUE)
# cnetplot(GSEA_enrichment, foldChange=mygene,node_label="category")#all,gene
# emapplot(GSEA_enrichment)



##2. Only_T_GSEA####
  sub_="Only_T"
  read_path <- paste0("./data/temp/",sub_,"_DEG_results_to_GSEA.csv")
  save_path_csv <- paste0("./data/temp/",sub_,"_GSEA_result_H.csv")
  save_path_polt <- paste0("./data/output/",sub_,"_GSEA_result_H_GSEA通路富集.pdf")
  read_path
  save_path_csv
  save_path_polt
  # 导入差异分析后的数据，以便后续使用logFC进行基因排序
  mydata <- read.csv(read_path ,row.names=1)
  head(mydata)
  mydata <- mydata[order(mydata$log2FoldChange,decreasing=TRUE),]#按log从大到小排序
  mydata <- mydata[mydata$log2FoldChange!=Inf & mydata$log2FoldChange!=-Inf ,]#去除无穷大的，不然后面报错，Not all stats values are finite numbers
  FCgenelist <- mydata$log2FoldChange
  names(FCgenelist) <- rownames(mydata)
  head(FCgenelist)
  # 开始GSEA富集分析
  GSEA_enrichment <- GSEA(FCgenelist,                 # 排序后的gene
                          TERM2GENE = geneSet_onco, # 基因集
                          pvalueCutoff = 0.05,      # P值阈值
                          minGSSize = 20,           # 最小基因数量
                          maxGSSize = 1000,         # 最大基因数量
                          eps = 0,                  # P值边界
                          pAdjustMethod = "BH")     # 校正P值的计算方法,可以做P矫正BH法，这样去掉假阳性
  #遇到报错重启RStudio解决：GSEA analysis... Error in serialize(data, node$con) : error writing to connection
  
  result <- data.frame(GSEA_enrichment)
  dim(GSEA_enrichment@result)
  write.csv(result, save_path_csv, row.names = FALSE)
  
  result$ID
  GSEA_enrichment@result$Description<- sub("HALLMARK_", "", GSEA_enrichment@result$ID)
  
  # 可视化
  # 特定通路作图——单个通路
  # p1 <- gseaplot2(GSEA_enrichment, "HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", color = "red3", pvalue_table = T) 
  # p1
  
  # 特定通路绘图——多个通路
  # p2 <- gseaplot2(GSEA_enrichment,
  #                 # result$ID,
  #                 c("HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION", 
  #                   "HALLMARK_TNFA_SIGNALING_VIA_NFKB",
  #                   "HALLMARK_HYPOXIA"),
  #                 #color = c("red3", "blue4","darkgreen"),
  #                 pvalue_table = T)
  # p2
  #单独通路的GSEA图
  # mygene<- GSEA_enrichment@geneSets[["HALLMARK_HYPOXIA"]]
  #   # c("LOX","ANGPTL4","SDC4","CCN1","COL5A1","GPC1","SDC2","BGN",
  #   #                  "GPC1",
  #   #                  "GPC3",
  #   #                  "GPC4")
  # p21 <- gseaNb(object = GSEA_enrichment,
  #        geneSetID = 'HALLMARK_HYPOXIA',
  #        newGsea = T,
  #        # addPoint = F,
  #        # newHtCol = c("blue","white", "red")
  #        addGene = mygene
  #        )
  # P21
  
  # mygene<- GSEA_enrichment@geneSets[["HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"]]
  # p22 <- gseaNb(object = GSEA_enrichment,
  #               geneSetID = 'HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION',
  #               newGsea = T,
  #               # addPoint = F,
  #               # newHtCol = c("blue","white", "red")
  #               addGene = mygene
  # )
  # p22
  # 展示富集到的通路，我们这里选择展示前15个
  p3 <- dotplot(GSEA_enrichment, 
                showCategory = 50,
                # color = "p.adjust",
                color = "p.adjust",
                # orderBy = 'NES',
                x = "NES",
                # x = 'enrichmentScore',
                title=sub_) 
  # p3
  
  # 将通路分为激活和抑制两个部分进行展示
  # 画图
  p4 <- dotplot(GSEA_enrichment, showCategory = 30, split = ".sign") + 
    facet_grid(~.sign) +
    theme(plot.title = element_text(size = 10, color = "black", hjust = 0.5),
          axis.title = element_text(size = 10,color = "black"), 
          axis.text = element_text(size = 10,color = "black"),
          axis.text.x = element_text(angle = 0, hjust = 1 ),
          legend.position = "right",
          legend.text = element_text(size = 10),
          legend.title = element_text(size = 10))+
    ggtitle(sub_)
  # p4
  
  #山脊图
  p5 <- ridgeplot(
    GSEA_enrichment,
    showCategory = 30,
    fill = "p.adjust",
    core_enrichment = TRUE,
    label_format = 30,
    orderBy = "NES",
    decreasing = FALSE)+
    ggtitle(sub_)
  # p5
  #气泡图-棒棒糖
  p6 <- dotplotGsea(data = GSEA_enrichment,
                    #topn = 30,
                    order.by = 'NES',
                    # line.type = 'dashed',
                    # pval=0.05, #defalut=0.05
                    # line.col = 'orange',
                    # add.seg = T
  )
  # p6
  p7 <- dotplotGsea(data = GSEA_enrichment,
                    #topn = 30,
                    order.by = 'NES',
                    # line.type = 'dashed',
                    # pval=0.05, #defalut=0.05
                    # line.col = 'orange',
                    add.seg = T)
  # p7
  
  save_path_polt <- paste0("./data/output/",sub_,"_GSEA_result_H_GSEA通路富集8x8.pdf")
  pdf(save_path_polt,width = 8,height = 8)
  # p1
  # p2
  print(p3)
  print(p4)
  print(p5)
  print(p6)
  print(p7)
  #不用printpdf打开显示文件损坏
  dev.off()
  save_path_polt <- paste0("./data/output/",sub_,"_GSEA_result_H_GSEA通路富集12x8.pdf")
  pdf(save_path_polt,width = 12,height = 8)
  # p1
  # p2
  print(p3)
  print(p4)
  print(p5)
  print(p6)
  print(p7)
  dev.off()

#barplot和dotplot都只显示最显著的富集项，而用户可能想知道哪些基因与这些显著项有关。为了考虑基因可能属于多个注释类别的潜在生物学复杂性，并提供可用的数值变化信息，我们开发了cnetplot函数来提取复杂的关联。cnetplot将基因和生物学概念(例如GO terms或KEGG pathways)之间的联系描述为一个网络。GSEA结果也只支持核心富集基因的表达。
# mygene<- GSEA_enrichment@geneSets[["HALLMARK_HYPOXIA"]]
# cnetplot(GSEA_enrichment, foldChange=FCgenelist, circular = TRUE, colorEdge = TRUE)
# cnetplot(GSEA_enrichment, foldChange=mygene,node_label="category")#all,gene
# emapplot(GSEA_enrichment)

